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We have performed a multicanonical molecular dynamics simulation on a simple model pro- 
tein. We have studied a model protein composed of charged, hydrophobic, and neutral spherical 
bead monomers. Since the hydrophobic interaction is considered to significantly affect protein 
folding, we particularly focus on the competition between effects of the Coulomb interaction 
and the hydrophobic interaction. We found that the transition which occurs upon decreasing 
the temperature is markedly affected by the change in both parameters and forms of the hy- 
drophobic potential function, and the transition changes from first order to second order, when 
the Coulomb interaction becomes weaker. 
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Protein structures and folding mechanisms are mysterious subjects^ According to Anfinsen's 
dogma,!' all of the information needed to construct a protein's three-dimensional structure is 
contained within its amino acid sequence. For given interactions between amino acids, we should 
in principle obtain a unique folding solution for transformation of proteins to their native states. 
However, from a numerical point of view, the prediction of a protein's native structure for a given 
amino acid sequence is considerably difficult because conformations accessible to a given polypeptide 
chain grow exponentially with chain length and it would require a much longer time than any 
computationally accessible time scale, known as the Levinthal (Levinthal paradox) 

In order to understand on a molecular level the general features of protein structures and mech- 
anism of folding, it is important to know the role of four major molecular interactions related to 
protein folding, which are hydrogen bonds, and van der Waals, hydrophobic, and Coulomb interac- 
tions. Computational studies for understanding protein folding mechanisms have been performed 
often under two particular assumptions, using either the simplest protein model with a lattice or 
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all-atom model together with solvent molecules in a continuous system. With these two extreme 
models, it is not easy to capture the nature of the intrinsic interaction in protein folding. To address 
this problem, one of the most promising methods is to examine simpler protein systems, such as 
bead-spring models, specifically designed to simulate the role of interactions. 

In general, complex systems with some different interactions often show frustration of potential 
energies, and thus exhibit a huge number of local-minimum-energy states in a free-energy landscape. 
Simulations of these systems by conventional molecular dynamics methods tend to get trapped 
in either of the local-minimum-energy states at low temperatures, and thus a long CPU time 
is needed to obtain accurate physical distributions since the relaxation time increases markedly 
with increasing complexity. Recently, numerical algorithms for overcoming such a multiple-minima 
problem have been proposed, such as the multicanonical ensemble method,^ which allows the 
trajectory of particles to escape over energy barriers and consequently enables sampling over a 
much wider phase space than the conventional method does. 

In this paper, we study a simple real protein, protein-g (PDB id: 2gbl), which is composed of 
56 amino acids. Figure 1 shows the distance map for protein-g, where filled symbols indicate the 
pairs of residues whose a-carbons are at a distance less than 13.0 [A]i Our model for protein-g is 
composed of positively (+) or negatively (— ) charged, hydrophobic (H), and neutral (P) spherical 
bead monomers with the following sequences: HPP + HHHPH + PH + H - PPP - HP - HHPH - 
+HH + PPHP - PHH - H - PPP - -HP + PHPHP- . Using this model for protein-g, we performed 
a multicanonical molecular dynamics (MMD) simulationjii) systematically. The main purpose 
of the present study is to examine the relative contributions among the respective interactions 
between monomers using molecular dynamics simulation, which may make it possible to construct 
the universal framework in protein folding. Since the hydrophobic interaction is considered to have 
a large effect on protein folding, it is of interest to observe how the resultant protein structure 
changes as the hydrophobic interaction sets in. Therefore, we employ the following interactions 
between monomers: (i) the excluded volume interaction (soft-core model), (ii) the covalent bond 
interaction (spring model), (iii) Coulomb interaction, and (iv) hydrophobic interaction. Recently, 
molecular dynamics simulations for the model with (i)-(iii) interactions have been performed by 
Baumketner et alS In the present study, the hydrophobic interaction, which originates from solvent 
effects, is also involved. Thus, in this paper we will in particular discuss the competition between 
effects of the Coulomb interaction and those of the hydrophobic interaction. Because of the lack of 
short-range attractive interactions such as hydrogen bond and van der Waals interaction, this model 
does not reproduce the true native structure of protein-g. However, since in this paper we mainly 
address conformational changes of a protein driven purely by long-range interactions, this is beyond 
the scope of the present study. Note that most theories of solution dynamics for homopolymers 
are based on the bead-spring model.!'© Advantages of the bead-spring model, compared to other 
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simplified models, such as lattice models, are that the model is practical, finds energy minima with 
less computation, and the physical properties obtained can be directly related to experiments. 

A solvent (like water) generally has the following two effects: (a) because of a screening effect, 
the interaction between charged monomers becomes significantly weaker (i.e., dielectric constant 
e* S> 1.0), (b) hydrophobic monomers repel water molecules in their neighborhood to result in a 
construction of hydrophobic cores and effectively gives rise to an attractive interaction (strength 
of the hydrophobic interaction is designated with e^b). Such solvent effects are non-trivial, because 
the number of associated solvent molecules is very large. In order to incorporate solvent effects 
into our model effectively, we tried to assume two different types of pair potential models for the 
hydrophobic interactions, the power potential V^ w ' and the exponential potential V^ p ' , which are 
given by 



where a v and d are 5.5 [A] and 10.0 [A], respectively. e sc is fixed at 3.0 [kcal/mol], and both a sc and 



In the case of the power potential, which was introduced by Shea et al.,ciP the interaction reaches 
only a few neighbors. On the other hand, in case of the exponential potential, which was used 
by Israelachvili et al., the attractive force reaches a much longer distance such as Coulomb 
interactions. As mentioned above, since the explicit form of the hydrophobic interaction is not yet 
well known, it is interesting to study the effects of the form of the hydrophobic interaction on the 
folding process. Therefore, we study here both types of potentials and compare the results obtained. 
In our simulation, we regard (e*,ehb) as the control parameters representing the strength of the 
Coulomb interaction relative to the hydrophobic interaction. A change of the control parameters 
corresponds to a change of solvents. We carried out various simulations including a two-limit 
system, i.e., pure vacuum and pure water (Table I). 

To analyze our data obtained by the present MD calculations, we used the following two order 
parameters, namely the radius of gyration (RG), which shows an order of spatial expansion of poly- 
mers quantitatively (eq. (§)), and the distance matrix error (DME), which measures the difference 
between structures of the model protein and a reference system rfj (eq. (||)); for the latter, we have 
chosen the native structure of protein-g: 
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Table I. The change of strength of the Coulomb interaction relative to the hydrophobic interaction in a two-limit 
system (i.e., pure vacuum and pure water). 



Interaction 


pure vacuum 


pure water 


Coulomb 


large (e* ~ 1) 


small (e* ~ 80.0) 


Hydrophobic 


small {tub ~ 0.0) 


large (ehb ~ 3.0[kcal/mol] ) 



DME 

where ./V is the number of residues, and r. 
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In MMD simulation, forces acting on particles are biased with a factor as shown below in eq. (5) 
which is calculated from the probability function of potential energy at each step so that we can 
obtain the flat distribution of potential energy P(E) after several preliminary simulations (eq. (^)). 



f(t) = -V£{E) 



d£(E) 
dE 



VT, 



-(1 + k B T — log P To (E))VE, 



£ ne W{E) = £ old {E) + kBTo l og p To ( E)i 



(5) 
(6) 



where k&, To, and Pt are the Boltzmann constant, a reference temperature and the canonical 
distribution function at temperature To, respectively. 

In Fig. 2, the probability distribution functions (PDF) of the potential energy during multicanon- 
ical runs and canonical runs (100, 200, and 300 [K]) are shown, in which the hydrophobic interaction 
obeys a power potential and the control parameters are (e*, ehb) = (5.0,3.0) at 300 [K]. The inset 
also shows the evolution of the potential energy during MMD simulation (10 million steps). Note 
that the time mesh in MD integration is fixed at 1 [fs] throughout this paper. It can be seen that 
the PDF through MMD after 30 iterations becomes flat over a wide range of the potential energy. 

Figure 3 shows the PDF obtained by MMD runs (and CMD runs in the inset) in two parameter 
planes, RG and DME, in which the hydrophobic interaction and the control parameter are the 
same as those in Fig. 2. It turns out that MMD runs samples more effectively than CMD runs do 
for the same number of samples (100 million)© 
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Using the reweighting methodHS to the probability distributions in MMD runs, we obtain canon- 
ical distributions for a wide range of temperatures. Figure 4 shows the temperature dependence of 
RG obtained in such a way. For the case that the hydrophobic interaction and the control parame- 
ter are the same as those in Fig. 2 (top-right panel of Fig. 4), we find two stable states with different 
RG values at approximately 300 [K] and a barrier between them, which suggest that the transition 
between these two states is of the first order. The corresponding distance maps for these two stable 
states are shown in Fig. 5. The right panel in Fig. 5 indicates that the C-terminal of our model 
protein-g is far from a folded state, in which RG takes a value as large as ~ 10 [A] at approximately 
300 [K]. When the hydrophobic interaction becomes weaker, for example (e*,e^) = (5.0,2.0) as 
shown in the top-left panel of Fig. 4, the transition seems to change from the first order to the 
second order. On the other hand, no transition exists for the case that the hydrophobic interaction 
is given by an exponential potential (bottom panels of Fig. 4). It follows that the existence of the 
transition changes markedly when the form of hydrophobic interaction is varied. 

To understand the influence of frustrations between several potentials more clearly, we employed 
weak Coulomb interactions to suppress frustrations. Figure 6 shows the temperature dependence 
of RG in weak Coulomb interactions, which corresponds to pure water (e*,^) = (80.0,3.0) as 
shown in Table 1. We found that the transition is of the second order instead of the first order at 
approximately 300 [K] despite the fact that the hydrophobic interaction has the same value as that 
in the top-right panel of Fig. 4. 

Here, let us consider the relationship between the present results and the experimental results 
for the gel phase transition. Polymers in the gel have only charged residues© For the gel phase 
transition to be of the first order, sufficient attractive and repulsive interactions are necessary. This 
suggests that the transition changes from the first order to the second order for the case that the 
Coulomb interaction becomes substantially weak. Our results for a model protein also show the 
same tendency for the first-order phase transition to be induced. In biopolymer systems, both 
Coulomb and hydrophobic interactions give rise to the energy frustration that plays a crucial role 
in the folding process. The character of transition is clearly changed when the solvent is varied© 
Our results suggest that the Coulomb interactions, both repulsive and attractive, contribute to the 
existence of the first-order transition. The general protein-folding mechanism is not yet known. 
Simulations to study the phase diagram for various model protein systems in various solvents are 
helpful and currently being undertaken by us. 

Finally we give below a brief summary of our work. We have performed multicanonical molecular 
dynamics simulations on a simple model for protein-g. Since we particularly focus on the competi- 
tion between effects of Coulomb and hydrophobic interactions, the control parameters (e*,ehb) are 
widely varied. We found that the present model for protein-g exhibits either continuous or discontin- 
uous phase transition, indicating explicit evidence of frustrations usually observed in heteropolymer 
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systems. RG dependence of temperature for each control parameter leads to the conclusion that 
the nature of the structural phase transition as a function of temperature is significantly affected 
by both the strength and the form of the hydrophobic interaction. 

The authors thank Dr. Andrij Baumketner for valuable discussions. This work was supported 
by the Research and Development for Applying Advanced Computational Science and Technology 
(ACT-JST). 
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Fig. 1. Distance map of protein-g (PDB id: 2gbl). The pairs of residues within 13.0 [A] are filled. The full sequence 
of the type of residues is also shown. 



6 



Fig. 2. The probability distribution functions of the potential energy for multicanonical runs (MMD) together with 
canonical runs (CMD) at 100, 200, and 300 [K]. The hydrophobic interaction obeys a power potential and the 
control parameters are (e*,e hb ) = (5.0,3.0) at 300 [K]. The inset also shows the evolution of potential energy 
during MMD simulation (10 million steps). 



Fig. 3. The probability distribution map in RG-DME plane obtained by multicanonical and canonical (inset) runs. 
The horizontal axis shows RG in the range between and 20 [A], and the vertical axis shows DME in the range 
between and 20 [A]. The hydrophobic parameters are the same as those in Fig. 2. 



Fig. 4. The temperature dependencies of the probability distribution of RG are shown for various control param- 
eters. The top figures are for the power hydrophobic interaction, and the bottom ones are for the exponential 
hydrophobic interaction. The horizontal axis shows RG in the range between and 20 [A], and the vertical axis 
shows temperature in the range of 200 ~ 600 [K] . 



Fig. 5. The distance maps for two stable states corresponding to the top-right panel of Fig. 4 at approximately 300 
[K]. 



Fig. 6. The temperature dependence of the probability distribution of RG in pure water is shown. The horizontal 
axis shows RG in the range of ~ 20 [A], and the vertical axis shows temperature in the range of 200 ~ 600 [K]. 
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